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Abstract 

A  simplified  dynamical  model  of  a  fuel  cell  of  the  proton  exchange  membrane  (PEM)  type,  based  on  physical-chemical  knowledge  of  the 
phenomena  occurring  inside  the  cell  has  been  developed  by  the  authors. 

The  model  has  been  implemented  in  the  MATLAB/SIMULINK  environment. 

Lab  tests  have  been  carried  out  at  ENEA’s  laboratories;  and  a  good  agreement  has  been  found  between  tests  and  simulations,  both  in  static 
and  dynamic  conditions. 

In  a  previous  study  [M.  Ceraolo,  R.  Giglioli,  C.  Miulli,  A.  Pozio,  in:  Proceedings  of  the  18th  International  Electric  Fuel  Cell  and  Hybrid 
Vehicle  Symposium  (EVS18),  Berlin,  20-24  October  2001,  p.  306]  the  basic  ideas  of  the  model,  as  well  as  its  experimental  validation  have 
been  published. 

In  the  present  paper,  the  full  implementation  of  the  model  is  reported  in  detail.  Moreover,  a  procedure  for  evaluating  all  the  needed 
numerical  parameters  is  presented. 

©  2002  Elsevier  Science  B.V.  All  rights  reserved. 
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1.  Introduction 

The  need  of  reducing  pollutant  emissions  and  of  utilising 
more  efficiently  the  available  energy  resources  (in  particular 
fossil  resources)  has  caused,  in  recent  years,  an  ever  increas¬ 
ing  attention  towards  fuel  cells.  In  fact,  their  high  conversion 
efficiency  and  low  environmental  impact,  make  them  good 
candidates  for  substituting,  at  least  in  some  applications, 
more  conventional  conversion  systems. 

One  of  the  applications  of  fuel  cells  currently  being 
considered  is  as  a  source  of  energy  for  electric  vehicles, 
normally  in  hybrid  configuration. 

Among  the  several  possible  kinds  of  fuel  cells,  the  polymer 
electrolyte  fuel  cells  (PEFCs),  also  called  proton  exchange 
membrane  fuel  cells  (PEMFCs),  appear  to  be  the  best  candi¬ 
date  for  the  use  aboard  of  electric  vehicles  in  which  simplicity, 
high  specific  power,  rapid  start-up  (even  at  low  temperature) 
have  the  maximum  importance  [4,11,16], 

To  be  able  to  utilise  these  devices  in  an  effective  way,  it  is, 
however,  mathematical  models  of  the  vehicle  fuel  cell  stack 
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are  necessary  so  that  the  system  behaviour  can  be  analysed  at 
the  design  stage  by  means  of  computer  simulations  in 
different  conditions  of  load,  pressure  of  reagent  gases, 
temperature. 

These  models  can  be  integrated  with  other  components  in 
vehicle  simulation  environments  (such  as  the  one  described 
in  [17]),  so  that  the  whole  system  of  a  fuel  cell  vehicle  can  be 
analysed  in  detail. 

Several  mathematical  model  of  PEM  fuel  cells  have 
already  been  presented  [4-12],  The  majority  of  them,  how¬ 
ever,  is  able  to  simulate  only  the  cell  steady-state  behaviour, 
while  the  analysis  of  their  performance  in  dynamic  condi¬ 
tions  is  important  for  the  use  aboard  of  vehicles,  given  the 
rapid  variation  of  mechanical  and  electrical  quantities. 

Some  models  [7,1 1]  are  characterised  by  a  high  complex¬ 
ity,  with  several  partial  differential  equations  to  be  kept  into 
account.  This  high  complexity  creates  problems  of  simula¬ 
tion  times,  parameter  identifications,  etc.  especially  when 
they  are  to  be  enclosed  into  a  larger  system,  such  as  the 
electric  vehicle. 

The  purpose  of  this  paper,  therefore,  is  to  propose  a 
dynamical  model  of  PEM  fuel  cells  that,  although  simplified 
and  containing  some  empirical  equations,  is  still  based  on 
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Nomenclature 

Ar  effective  Pt  surface  area  per  unit  geometric 

surface  area  (cm2/cm2) 

Am  total  Pt  surface  area  per  unit  geometric  surface 
area  (cm2/cm2) 
b  Tafel  slope  (V) 

cH+  dimensionless  proton  concentration 

ck  concentration  of  species  k  (mol/cm3) 

Qi  double-layer  capacitance  (F/cm2) 

3sik  total  pressure-diffusivity  product  (bar  cm2/s) 

Dik  diffusivity  of  the  gas  pair  i-k  (cm2/s) 

Ea  open-circuit  voltage  (V) 

EoC  cathode  potential  (V) 

£Vef  reference  voltage  (V) 

/  frequency  (s_1) 

F  Faraday’s  constant  (96484  C/eq.) 

[H+]  proton  concentration  (mol/cm3) 

j  cell  current  density  (A/cm2) 

j0  exchange  current  density  (A/cm2) 

jT  reaction  density  current  (A/cm2) 

Ld  thickness  of  diffusion  layer  (cm) 

/Vair  inlet  air  mole  flow  (mol/(cm2  s)) 

Nh2  hydrogen  consumption  (mol/(cm2  s)) 

Nk  superficial  flux  of  species  k  (mo]/(ctrr  s)) 
Aide  oxygen  flux  at  interface  diffusion  layer- 

catalyst  layer  (mol/(cm2  s)) 
pc  cathode  total  pressure  (atm) 

pk  partial  pressure  of  species  k  (bar) 

Pic  oxygen  partial  pressure  within  catalyst  layer 
(bar) 

Pio  zero-current  oxygen  partial  pressure  (bar) 

R  gas  constant  (82.056  (bar  cm3)/(K  mol)  = 

8.3144  J/(K  mol)) 

Rohm  internal  ohmic  resistance  (ohm  cm") 

T  absolute  temperature  (K) 

xk  mole  fraction  of  species  k 

y  distance  through  diffusion  layer  (cm) 

Greek  letters 

fl  inlet  nitrogen-oxygen  mole  ratio 

Eg  porosity  of  diffusion  layer 

A(/jce  catalyst-electrolyte  potential  difference  (V) 

r\  cathode  overpotential  (V) 

t  tortuosity  of  diffusion  layer 

th+  time  constant  of  proton  concentration  (s) 

£  dimensionless  distance 


the  chemical-physical  knowledge  of  the  phenomena  occur¬ 
ring  inside  the  cell. 

The  effort  carried-out  was  to  simplify  as  much  as  possible 
the  analytical  aspects,  so  that  simulation  and  numerical 
parameter  identification  are  eased. 

The  model  is  implemented  in  the  well  known  and  highly 
widespread  simulation  environment  MATLAB/SIMULINK. 


2.  Mathematical  definition  of  the  proposed  model 

2.1.  A  simple  description  of  the  cell  structure 

A  PEM  fuel  cell  is  constituted  by  a  membrane  able  to 
conduct  protons  disposed  between  two  electrodes.  The 
electrodes-membrane  assembly,  in  turn,  is  pressed  by  two 
conductive  plates  containing  some  channels  in  which  the 
reactants  flow. 

A  simplified  representation  of  the  cell  is  reported  in  Fig.  1 . 
The  main  elements  composing  the  cell  are:  conductor  plates, 
electrodes  and  membrane. 

The  electrodes  are  constituted  by  a  gas  diffusion  layer 
(made  of  carbon  paper  or  cloth,  carbon  powder  and  PTFE) 
and  a  catalyst  layer  (constituted  by  Pt/C  powder  and  Nafion 
ionomer);  both  have  a  porous,  partially  hydrophobic,  structure. 

According  to  several  studies  [3-11],  the  larger  pores 
(called  macropores)  operate  as  a  ducts  for  the  reactants 
from  the  flow  channels  towards  the  catalyst  layer,  while 
the  smaller  ones  ( micropores )  operate  as  a  ducts  for  the 
passage  of  water. 

The  reactions  occurring  in  the  catalytic  layer  of  the 
electrodes  are  the  following  ones:  at  the  anode  the  hydrogen 
decomposes  and  yields  electrons  to  platinum  and  protons  to 
the  membrane: 

H2  -f  2H+  +  2e“  (1) 

at  the  cathode,  the  oxygen  reacts  with  the  protons  coming 
from  the  membrane  and  with  the  electrons  supplied  by  the 
catalyst  and  forms  molecules  of  water: 

5O1  +  2H+  +  2e“  — >  H20.  (2) 

The  membrane,  if  well  humidified,  has  a  proton  conductivity 
sufficiently  high. 

Inside  the  cells,  however,  because  of  the  dragging  of  water 
molecules  by  the  proton  flow,  the  part  of  the  membrane  on 
the  cathode  side  tends  to  saturate,  while  the  one  on  the  anode 
side  tend  to  dehydrate,  with  a  consequent  conductivity 
reduction;  to  limit  this  phenomenon,  and  to  avoid  water 
loss  by  evaporation,  the  cell  is  fed  with  humidified  gases  at 
higher  temperatures  than  the  cell  temperature. 

2.2.  The  cell  voltage  versus  current  density  behaviour 

The  PEM  fuel  cell  behaviour  is  quite  complex  and  is 
influenced  by  several  factors,  including  composition  and 
structure  of  electrodes  and  membrane  type. 

An  analysis  of  a  large  amount  of  experimental  data 
reported  in  literature  [7,11,14,23-25]  has  shown  that  the 
voltage  versus  current  density  plots  for  a  PEM  cell  have  a 
behaviour  of  the  type  shown  in  Fig.  2,  in  which  the  numer¬ 
ical  values  are  quite  typical  although,  obviously,  the  actual 
values  of  any  particular  cell  can  differ  from  them. 

Since  the  early  1960s,  a  lot  of  papers  have  been  published 
on  the  modelling  of  PEM  cells.  The  objectives  of  these 
studies  were  to  analyse  the  effects  of  composition  and 
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Fig.  2.  Typical  cell  voltage  vs.  current  density  plots  for  PEM  fuel  cells  and  a  common  interpretation  for  the  voltage  drop. 


structure  of  cell  on  the  cell  or  half-cell  voltage  versus  current 
density  behaviour  [9]. 

Initially,  less  attention  was  given  to  cathode  problems  than 
was  devoted  to  ionomeric  membranes.  This  initial  emphasis 
on  voltage  drops  in  the  membrane  electrolyte  was  required 
because  this  unique  component  of  cell  is  quite  different  from 
the  electrolytes  employed  in  other  low-temperature  fuel  cell 
types  [7],  With  the  development  of  membrane  and  particularly 
with  the  introduction  of  thinner  types,  it  has  become  clear  that 
voltage  drops  in  membrane  of  a  well-humidified  PEM  cell  are 
not  always  the  prevalent  ones,  particularly  in  cells  with  air  as 
oxidant.  With  membranes  as  Nafion  112  (50  pm  thickness), 
Membrane  Dow  (80  pm  thickness),  Membrane  “C”  (140  pm 
thickness),  the  cell  resistance  in  well-humidified  PEM  cells 
does  not  increase  excessively  with  the  current,  even  at  den¬ 
sities  as  high  as  2.5-3  A/cm2  [7,14], 

According  to  the  already  mentioned  literature,  the  voltage 
drops  in  a  PEM  fuel  cell  can  originate  from 

1.  limited  cathode  interfacial  kinetics,  which  determines 
high  voltage  drops  mainly  in  low  current  density  region 
(Fig.  2),  both  with  air  and  oxygen  cathode; 


limited  protonic  conductivity  in  the  catalyst  layer  and  in 
the  membrane  electrolyte;  in  the  case  of  oxygen  cathode 
the  ohmic  drops  within  the  membrane  constitute  the  major 
fraction  of  measured  drops  in  the  current  range  of  interest; 
anode  interfacial  kinetics,  which  determines  a  voltage 
drop  linear  with  the  current  density  because  of  the  high 
exchange  current  density  for  the  electrooxidation  of 
hydrogen;  however,  these  drops  result  to  be  very  small 
when  compared  to  the  others; 

limited  cathode  mass  transport,  particularly  in  the  case 
of  air  cathode,  which  determines  a  rapid  decrease  of  cell 
voltage  with  current  density;  this  diffusion  limitation 
lead  to  pseudo  resistive  effects  at  moderate  current 
densities  and  the  appearance  of  limiting  currents  at  high 
current  densities; 

cathode“flooding” ,  when  the  water  produced  from 
cathode  reaction  is  not  removed  effectively  from  the 
cell;  consequently  the  excess  liquid  water  leads  to  a  poor 
transport  of  oxygen  through  the  electrode;  this  phenom¬ 
enon,  with  optimised  cells,  occurs  typically  in  oxygen 
cathode  at  very  high  current  densities,  which  cannot  be 
reached  if  the  cell  cathode  is  fed  with  air. 
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2.3.  Basic  assumptions  and  model  structure 

The  model  considers  a  cell  that  utilises  H2  as  a  fuel  and  air 
as  oxidant,  both  humidified. 

Based  on  the  basic  considerations  reported  in  the  previous 
section,  the  main  assumptions  of  the  proposed  model  are  as 
follows: 

1.  The  electrode  (diffusion  and  catalyst)  layers  have  been 
considered  as  having  given,  and  constant  porosity  and 
tortuosity.  This  implies  that  the  possibility  of  cathode 
flooding  is  not  considered,  since  it  would  change  the  size 
of  macropores  and  therefore  the  layer  effective  porosity. 

2.  The  model  is  one-dimensional,  i.e.  all  quantities  vary 
only  in  the  direction  orthogonal  to  anode  and  cathode 
surfaces. 

3.  The  temperature  is  supposed  to  be  uniform  in  the  cell. 

4.  As  often  made  by  authors  [5-10]  the  air  total  pressure  is 
assumed  to  be  uniform,  while  the  variation  with  space  of 
the  partial  pressures  of  its  components  is  kept  into 
account.  As  a  consequence  of  this,  the  movement  of 
gases  within  electrodes  is  consequence  of  just  concen¬ 
tration  and  not  pressure  gradients. 

5.  The  water  vapour  contained  in  the  reactants  of  the 
baking  macropores  is  in  equilibrium  with  the  surround¬ 
ing  liquid  phase;  consequently  the  partial  water  pressure 
is  uniform. 

6.  The  membrane  is  considered  as  being  completely 
saturated  of  water;  therefore,  its  conductivity  is  only 
function  of  temperature. 

7.  The  anodic  overpotential  is  disregarded  with  respect  to 
the  cathode  ones  (cf.  [7,8,16]);  consequently,  the  voltage 
across  the  anode  is  considered  to  be  constant. 

According  to  these  assumptions,  the  proposed  mathema¬ 
tical  model  is  composed  by  the  elements  shown  in  Fig.  3,  in 


which  also  the  logical  connections  of  the  different  blocks, 
and  the  model  inputs  and  outputs  are  shown. 

In  the  following  paragraphs,  a  description  of  the  imple¬ 
mentation  of  the  different  blocks  of  Fig.  3  is  reported  making 
reference,  as  usual,  to  current  densities  instead  of  actual  cell 
and  reaction  current. 

2.4.  Cathode  gas  diffusion 

The  oxygen  contained  in  the  air  entering  the  cell  before 
reaching  the  catalyst  layer  diffuses,  through  the  diffusion 
layer,  within  a  gaseous  mix  constituted  by  nitrogen,  water 
vapour  and  the  oxygen  itself,  all  considered  as  being  ideal 
gases.  This  phenomenon  is  described  [5,8]  by  the  following 
differential  equations: 

continuity  equations: 


h_d_Pi  dNi  = 

RT  dt  dy 

Stefan— Maxwell  equations  of  diffusion: 


Eg  d pi 

x2  dy 


E 


iPc^ik 


(. PiNk-PkNi ), 


(3) 


(4) 


where  i,  k  £  (1,  3)  and  p\  is  the  oxygen  partial  pressure, 
Pi  =  Psat  ( T)  the  water  vapour  partial  pressure,  p3  the  nitro¬ 
gen  partial  pressure  and  pc@ik  =  Oik  =  DiffT),  in  which 
pc=pi+p2+  p3  is  the  cathode  (air)  pressure,  assumed,  in 
accordance  to  hypothesis  4  stated  in  Section  2.3,  constant. 

The  analytical  expressions  of  the  diffusivities  Dik{T)  have 
been  taken  from  [13]. 

The  diffusion  layer  porosity  6g  (ratio  between  pore  volume 
and  total  layer  volume)  is  considered  independent  on  the 
cell  operating  conditions;  the  parameter  x  (tortuosity:  ratio 


Fig.  3.  The  general  model  structure. 
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between  actual  pore  length  and  macroscopic  diffusion 
layer  thickness)  keeps  into  account  the  fact  that  the  distance 
the  oxygen  covers  for  reaching  the  catalytic  layer  is  greater 
than  the  diffusion  layer  thickness,  because  of  the  pore 
curvature. 

Combining  Eqs.  (3)  and  (4)  and  applying  the  hypotheses 
stated  in  the  previous  paragraph,  the  following  single  partial 
differential  equation  can  be  obtained,  having  the  unique 
unknown  p\  (oxygen  partial  pressure): 


dpi  _  d2pi  jv  dpi 

dt~md £  ^4 


(5) 


in  which  q  is  the  dimensionless  abscissa  ylLd,  the  variable 
/Vidc  is  the  oxygen  flow  at  the  interface  diffusion-catalyst 
layers,  and  a>  and  t ji  are  function  of  cell  temperature  and  air 
pressure  only: 


^Lliipwt/Dn)  +  (Pc-Ps^/Dn)  ’ 


<A 


RT 

SghdOc  -  Psat)  ‘ 


(6) 


In  Appendix  A,  mathematical  details  of  the  analysis  by 
which  Eq.  (5)  is  derived  are  presented. 


2.5.  Cathode  kinetics 


The  electrochemical  reaction  that  occurs  in  the  catalyst 
layer  is  the  (2),  in  which  the  oxygen  reacts  with  the  protons 
present  in  the  electrolyte  and  with  the  electrons  supplied  by 
the  catalyst,  thus  forming  water  molecules. 

The  current  that  is  created  by  the  electrochemical  reac¬ 
tion,  in  addition  to  the  reactants  concentration,  depends  on 
the  potential  difference  between  catalyst  and  electrolyte,  as 
described  by  the  Butler- Volmer  equation  [9,11]: 


t=W^Kexp©“1} 

1  =  EoC  -  Ar/jcc  =E0-  Vceii  -  Rdhmj. 


(7) 

(8) 


The  cell  current  density,  in  turn,  is  the  sum  of  the  reaction 
current  plus  the  contribution  due  to  the  charge  storage  in  the 
electrical  double-layer  at  the  interface  between  catalyst  and 
diffusion  layers: 


j  —  jt  +  Qi  v 


the  double-layer  capacitance  Cdl  is  assumed  to  be  constant 
[16]. 

Some  previous  studies  [5,6,11,12]  have  shown  that  the 
quantity  of  platinum  interested  by  the  chemical  reaction 
lowers  as  the  cell  current  rises;  i.e.  the  reaction  tends  to 
concentrate  in  a  smaller  region  of  the  catalyst  layer,  parti¬ 
cularly  when  the  cell  is  fed  with  air.  This  has  been  attributed 
to  the  limited  electrolyte  conductivity  and  limited  oxygen 


transport  that  determine  voltage  drop  and  02  concentration 
gradient  respectively. 

To  keep  into  account  this  phenomenon,  without  introdu¬ 
cing  a  large  number  of  partial-derivative  differential  equa¬ 
tions,  the  parameter  Ar  (ratio  between  total  catalytic  surface 
and  cell  surface)  has  been  considered  to  be  function  of  the 
cell  current  density,  according  to  the  following  equation: 

Ar  =  Aro  exp(— aij  -  ay5)  (10) 

in  which 

ai  „,2  (h  ki  \ 
ai  =  ai.K(Po2N^  exp  I  —  -  —  I , 

Psat.,a  3  fk2  k2\  .... 

a2  =  «2,ref-^exp^---j.  (11) 

The  proposed  expressions  of  a\  and  a2  have  been  defined  to 
interpolate  at  best  the  experimental  results  obtained  at 
different  temperatures  and  pressures  of  the  reactant  gases 
(see  Section  3.4).  Eqs.  (10)  and  (11)  are  an  empirical 
mathematical  way  to  overcome  the  high  computational 
complexity  of  the  partial  differential  equations  required  to 
model  the  above-mentioned  phenomenon. 

In  addition  to  (10),  the  variables  that  appear  in  (7)  are 
considered  to  be  functions  of  time  only. 

From  experimental  measures  related  to  the  reduction  of 
oxygen  in  an  acid  environment  [1,2]  it  is  known  that 
parameters  b  and  ja  of  (7)  present  a  sudden  rise  for  values 
of  cathode  voltage  of  about  0.80  V. 

In  [1],  quantitative  data  on  temperature  dependence  of 
electrode  kinetics  of  oxygen  reduction  at  platinum/Nafion 
interface  are  reported.  All  of  the  Tafel  plots  reported  show 
two  distinct  slopes  at  all  temperatures  investigated  and  the 
break  in  the  plots  occurs  at  a  potential  of  0.75-0.8  V.  This 
break  has  been  observed  previously  and  interpreted  as  a 
change  in  oxygen  reduction  mechanism  [2],  The  low 
Tafel  slope,  corresponding  to  a  potential  regime  where 
oxygen  reduction  proceeds  on  Pt-oxide  covered  surface, 
results  dependent  on  temperature  according  to  the  following 
relation 


where  a  =  1  and  n=  1  has  been  assumed. 

The  high  Tafel  slope,  corresponding  to  a  potential  regime 
where  oxygen  reduction  proceeds  on  free-oxide  Platinum, 
does  not  show  a  monotone  behaviour  with  the  temperature 
and  can  be  considered  temperature  independent;  from 
experimental  data  reported  in  [1]  with  the  temperature 
changing  from  30  to  80  °C  the  following  constant  value 
can  be  obtained: 


bh  =  0.0504  V.  (13) 

The  exchange  current  density  ja  shows  a  dependence  both  on 
temperature  and  oxygen  partial  pressure,  and  as  the  Tafel 
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slope,  at  high  current  densities  it  assumes  different  values 
from  ones  at  low  current  densities. 

To  keep  into  account  in  the  model  of  above-mentioned 
behaviour  of  the  Tafel  slope  and  the  exchange  current 
density,  the  two  parameters  ( b  and  ja)  has  been  defined  as 
a  function  of  the  cathode  overpotential  r y. 

jo  =  [/oi  +  A /  •  u(t]  -  rib)]  (14) 

b  =  h+{b h  -  bx)  ■  u(r,  -  r/b )  (15) 

in  which  m(  )  is  the  Heaviside  function  and 
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Fig.  4.  Cell  dynamic  response  to  a  current  step. 


r,b  =  E0-  0.8  (16) 

4/=J„,[exp(|-^)-l]  (17) 

the  overpotential  rjb  is  the  cathode  overpotential  correspond¬ 
ing  to  break  point  in  the  Tafel  plots  on  oxygen  reduction  at 
platinum/Nafion  interface. 

The  A /'  expression  in  (17)  has  been  obtained  imposing  the 
continuity  of  (7)  with  r]  =  rjb. 

Considering  Eqs.  (10),  (14)  and  (17),  the  product  joAT  that 
appears  in  (7)  becomes: 


joAr 


x  exp (-cnj-aij5)  (18) 


in  which,  as  shown  also  in  literature  [5,6],  the  product  j„\Aro 
results  to  be  function  of  temperature  and  oxygen  partial 
pressure: 

JolAro  =  J  ref  exp  j  +  f^Po2  ■  (19) 

Differently  from  what  assumed  in  other  papers  [4,6,9,11], 
the  proton  concentration  [H+]  in  the  catalyst  layer  has  been 
considered  to  be  a  function  of  the  cell  current.  This  because 
if  the  cell  current  rises  the  water  production  and  the  hydra¬ 
tion  of  the  polymeric  electrolyte  near  the  cathode  catalyst 
tend  to  rise  as  well,  causing  this  way  an  increase  in  the 
number  of  mobile  protons.  In  fact,  from  experimental  data 
[15],  it  is  known  that  the  proton  conductivity  of  a  polymer 
electrolyte  increases  with  the  water  content  in  the  polymer, 
and  this  behaviour  is  attributed  to  the  increase  of  the  number 
if  mobile  protons  with  the  water  content. 

For  a  polymer  electrolyte  the  conductivity  is  defined  as 
[5,10]: 


a 


F2 

RT 


Dn+cH+ 


(20) 


so  if  it  increases  with  the  water  content  the  proton  concen¬ 
tration  increases  as  well,  according  to  (20). 

This  hypothesis  is  able  to  explain  the  voltage  overshoot 
visible  in  the  cell  experimental  response  to  sudden  current 
reductions  (Fig.  4),  that  was  already  observed  in  [18]  with¬ 
out  explanations:  if  the  cell  current  increases  very  fast  the 
proton  concentration  near  the  cathode  catalyst  increases 


very  fast  as  well,  while  if  the  cell  current  decreases  very 
fast  the  proton  concentration  decreases  more  slowly  follow¬ 
ing  the  slow  dynamics  of  the  excess  water  removal  from  the 
cathode  catalyst  layer. 

Therefore,  the  behaviour  over  time  of  the  proton  concen¬ 
tration  [H+]  has  been  defined  by  the  following  empirical 
differential  equation: 


(  gfgA  dcH+  cH+  _  1  +«H+  •/ 

V  dt  )  dt  th+  %+ 


(21) 


in  which  u(  )  is  the  Heaviside  function,  CH+  =  [H+]/[H+]0 
the  dimensionless  proton  concentration,  th+  the  time  con¬ 
stant  related  to  the  CH+  dynamics,  and  aH+  is  a  parameter 
that  links  CH+  to  the  cell  current. 


2.6.  Internal  resistance  and  voltage  drops 


The  passage  of  current  through  the  cell  causes  ohmic 
voltage  drops  basically  due  to  the  electron  transfer  in 
electrodes  and  in  the  conductive  graphite  plates  and  to 
the  proton  transfer  through  the  membrane. 

Since  the  conductivity  of  graphite  is  much  larger  than  the 
one  of  the  membrane  [4,12],  the  drops  due  to  the  electron 
transfer  can  be  neglected;  therefore,  the  cell  internal  resis¬ 
tance  practically  coincides  with  that  of  membrane. 

In  general,  the  membrane  resistance  results  to  be  a  func¬ 
tion,  in  addition  of  temperature,  also  of  the  current  [14,15], 
As  shown  in  some  previous  studies  related  to  the  simulation 
of  the  Nafion  membrane  behaviour  [20,21],  when  the  current 
rises  the  dragging  from  anode  to  cathode  of  water  molecules 
by  the  current  causes  a  rising  concentration  gradient.  As  a 
consequence,  the  side  of  the  membrane  on  the  cathode  and 
anode  parts  tends  to  saturate  and  dehydrate  respectively, 
with  a  reduction  of  the  local  conductivity  and  an  increase  in 
the  membrane  resistance.  This  phenomenon,  with  optimised 
cells,  occurs  typically  in  oxygen  cathode  at  very  high  current 
densities,  which  cannot  be  reached  if  the  cell  is  fed  with  air. 

As  far  as  the  influence  of  temperature  on  membrane 
conductivity,  it  can  be  said  that,  for  temperatures  below 
100  °C,  if  the  temperature  increases  the  conductivity 
increases  as  well. 

In  [15],  experimental  data  on  the  temperature  dependence 
of  proton  conductivity  of  some  membranes  are  reported.  The 
tested  membranes  were  Nafion  117,  “C”  Membrane  and 
Dow  Membrane,  and  their  conductivity  was  determined 
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Fig.  7.  Internal  resistance  at  60  °C  and  different  pressures  of  reactant  gas. 


both  with  membranes  immersed  in  water  and  with  mem¬ 
branes  partially  hydrated.  All  membranes  immersed  in  water 
showed  a  quite  linear  behaviour  of  the  conductivity  with  the 
temperature  over  the  range  30-90  °C.  In  Fig.  5,  conductivity 
values  of  Nafion  117  are  reported. 

In  [22],  the  proton  conductivity  of  Nafion  117  was 
measured  under  various  conditions  of  humidity  and  tem¬ 
perature,  and  all  measures  were  carried  out  with  the  mem¬ 
brane  under  exposure  to  a  stream  of  humid  gas.  The 
conductivity  measured  at  different  temperatures  exhibited 
a  behaviour  as  shown  in  Fig.  6.  In  this  case,  the  conductivity 
has  not  a  monotone  behaviour  with  the  temperature,  but 
from  20  to  45  °C  it  decreases  and  then  it  increases  when  the 
temperature  increases  from  45  to  80  °C.  This  behaviour  was 
attributed  to  the  decrease  of  water  content  in  the  membrane 
between  20  and  45  °C,  considering  that  above  45  °C  the 
water  content  was  rather  constant. 

From  the  experimental  observations  made  by  the  ENEA’s 
labs  related  to  a  PEM  cell  having  a  Nafion  115  membrane, 
the  internal  resistance  resulted  to  be,  in  practice,  function  of 
only  temperature  in  the  range  from  40  to  80  °C.  In  Figs.  7 
and  8,  some  values  of  cell  resistance  at  60  and  70  °C  are 
reported.  The  resistance  was  determined  as  high  frequency 
value  of  cell  impedance  for  small  ac  signal  around  steady- 
state  current  values  (Section  3.2).  As  shown  in  the  figures  the 
resistance  resulted  to  be  independent  on  gas  pressures; 
moreover,  as  the  temperature  changes,  the  behaviour  of 
resistance  versus  current  may  switch  from  decreasing  to 
increasing. 
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Fig.  6.  Proton  conductivity  of  Nafion  117  under  exposure  to  a  stream  of 
humid  gas  (from  [22]). 


As  a  consequence,  for  sake  of  simplicity,  at  each  tem¬ 
perature  the  cell  resistance  was  chosen  as  the  average  of  all 
values  obtained  at  the  different  currents  and  gas  pressures. 

In  Fig.  9,  the  resistance  values  determined  in  the  range 
from  40  to  80  °C  are  reported;  in  the  model  they  have  been 
interpolated  by  the  following  equation: 

Roim  =  RKf  +  *T{T-TKi)  (22) 

in  which  Trc(  =  343.15  K. 

In  Fig.  10,  a  comparison  of  conductivity  values  deter¬ 
mined  in  this  work  with  those  reported  in  [22],  and  obtained 
with  the  membrane  operating  in  conditions  similar  to  those 
occurring  in  the  actual  cell,  is  shown. 

The  comparison  shows  a  good  agreement  in  terms  of  trend 
of  conductivity  with  the  temperature;  the  marked  differences 
in  terms  of  numerical  values  shown  in  Fig.  10  have  little 
significance  since  important  differences  exist  in  the  measure 
conditions  of  the  two  cases. 
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Fig.  10.  Comparison  between  the  conductivities  utilised  in  this  work  and 
in  [22], 


This  result  can  be  obtained  considering  that,  in  steady-state 
conditions,  ja  and  b  can  be  expressed  as  a  function  of  j: 
jo  =  joi  +  4/  •  u(j  -jb),  b  =  b\  +  (bh  —  b{)  ■  u(j  —  jb ) 

(28) 

in  which  jb  is  the  current  density  that,  in  steady-state, 
corresponds  to  r)b. 

The  Eq.  (27)  is  used  in  this  work  to  determine  the 
parameters  j0|,  Am  ax,  a2  and  v  as  functions  of  cell  tem¬ 
perature,  gas  pressures  and  gas  flows  (Section  3.4). 


2.7.  Open-circuit  voltage  (OCV) 

The  open-circuit  voltage  has  been  defined  by  means  of  the 
Nernst  equation 

Eo  =  Eref  +  ^  (T  -  Tier)  +  g  HpbJ>£)  (23) 

applied  to  cell  reaction 

H2(g)+i02(g)^H20(aq). 

To  improve  agreement  to  the  experimental  results,  the 
Eq.  (23)  is  slightly  modified  as  follows: 

Eo  =  Eref  +  ^(E-  Tief)  +  \n(pHlp{£)  (24) 

in  which  the  parameter  k  has  been  introduced.  In  fact,  with  a 
k  of  1  and  standard  values  of  Eref,  dEldT,  the  E0  would  be 
much  higher  than  the  values  experimentally  observed,  as 
already  noted  in  the  past  [1], 

The  reference  temperature  Tref  has  been  fixed  to  343.15  K. 


3.  Identification  of  numerical  values  of  parameters 

The  proposed  model  contains  22  numerical  parameters 
that  need  to  be  identified;  they  are  shown  in  Table  1,  along 
with  the  equations  in  which  they  appear. 

Among  which,  the  parameters  ft  and  bh  can  be  directly 
found  in  literature  ([5]  and  [1],  respectively).  In  the  follow¬ 
ing  a  possible  procedure  for  determining  all  the  other 
parameters  of  Table  1  is  presented. 

3.1.  Open-circuit  voltage  (parameters  Err.j,  dEJdT,  k) 

The  cell  voltage  stabilises  completely  after  several  sec¬ 
onds  (a  minute  or  so)  from  any  change  in  the  parameters 
from  the  exterior,  i.e.  reactant  pressures,  cell  temperature 
and/or  current  density. 

If  the  OCV  is  measured  for  different  values  of  pn2  ,Po2  and 
T,  the  three  parameters  Eief,  AEJAT,  k  can  be  easily 
found  using  Eq.  (24).  In  Table  2,  some  experimental  values 
of  OCV  at  different  temperatures  and  pressures  are  reported. 


2.8.  Steady-state  cell  voltage 

The  simplifying  assumption  contained  in  Eq.  (10)  allows 
to  derive  an  analytical  expression  able  to  fit  the  experimental 
data  of  cell  voltage  over  the  whole  range  of  current. 

The  steady-state  solution  of  Eq.  (5)  (Appendix  A)  gives 
the  oxygen  partial  pressure  in  the  diffusion  layer  as  function 
of  cell  density  current: 

Ei  (£)  =  Pio  [1  +  P  —  P  exp(y/£)]  (25) 

in  which  v  =  \j//AtoF  then  the  oxygen  partial  pressure  in  the 
catalyst  layer  (<?  =  1)  will  be 

Ei(l)  =Eio[l  +  /?-£exp(y/)].  (26) 

Combining  Eqs.  (7)-(9)  and  (26)  the  following  analytical 
expression  of  the  cell  voltage  can  be  obtained: 

Vceii  =  Eo  +  b  ln(/cAo)  -  (axj  +  azj5) 

+  Mn|[1  +  g  — ^,pM,(1+w-3,|_ 

l  J  +JoAr  J 

(27) 


Table  1 


Numerical  parameters  of  the  model 


Parameter 

Equation 

rid 

(6) 

6gid 

(6) 

P 

(25) 

Cdi 

(9) 

aijef 

(11) 

ki 

(11) 

(11) 

Oil 

(11) 

<J2jef 

(11) 

k2 

(11) 

a3 

(11) 

bh 

(15) 

/ref 

(19) 

ka 

(19) 

«0 

(19) 

th+ 

(21) 

“h+ 

(21) 

Aef 

(22) 

aT 

(24) 

£rc, 

(24) 

dEJdT 

(24) 

k 

(24) 
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Table  2 

Variation  of  open-circuit  voltage  with  temperature  and  pressures  of 
reactants 


T(°C)  pc  (bar) 

Pa  (bar) 

Po2  (bar) 

Ph2  (bar) 

OCV  (V) 

70  1 

1 

0.143 

0.680 

0.959 

70  2 

1 

0.349 

0.680 

0.966 

70  3 

2 

0.556 

1.667 

0.978 

70  4 

3 

0.763 

2.654 

0.984 

50  3 

2 

0.594 

1.851 

0.985 

30  3 

2 

0.611 

1.930 

0.989 

Table  3 

Numerical  parameters  of  model 

Parameter 

Unit 

Value 

rLd 

cm 

0.149 

£gZ/d 

cm 

0.026 

P 

- 

3.773 

Cat 

mF/cm2 

41.71 

alref 

cm2/A 

25,460 

ki 

K 

1174 

at 

- 

0.238 

a2 

- 

0.805 

d2,ref 

(cm2/A)5 

0.127 

k2 

K 

20160 

«3 

- 

-0.566 

K 

V 

0.0504 

jre  f 

A/cm2 

3.875  x 

:  10  3 

ka 

K 

18300 

aD 

- 

0.583 

th+ 

s 

12.78 

aH+ 

(cm2/A)3 

5.87 

/?ref 

flcm2 

0.269 

aT 

(SI  cm2)/K 

2.02  x 

10  3 

#ref 

V 

0.975 

dEJdT 

mV/K 

0.27 

k 

- 

0.755 

The  partial  pressures  of  02  and  H2  have  been  determined 
considering  at  each  temperature  the  water  partial  pressure. 

The  interpolation  of  experimental  data  by  Eq.  (24)  yields 
the  numerical  values  of  parameters  Ere  f,  dEJdT,  k  reported  in 
Table  3. 

The  Fig.  1 1  shows  the  comparison  between  experimental 
and  calculated  values  after  the  fitting. 


-0.5  0  0.5 

^{pHlPo]) 


Fig.  11.  Experimental  and  model  values  of  open-circuit  cell  voltage  at 
70°C  and  different  reactant  pressures. 


Fig.  12.  (a)  AC  impedance  spectrum  at  low  overpotentials;  (b)  Equivalent 
circuit  of  ac  impedance  spectrum  at  low  overpotentials. 


3.2.  Internal  impedance  ( parameters  Rrej,  txT,  Cdi) 


If  the  cell  is  stimulated  with  small-amplitude  signals  at 
various  frequencies  a  diagram  of  the  complex  impedance 
can  be  drawn  in  the  Gauss  plane. 

Experimental  results  reported  in  [16]  show  that  this 
impedance,  while  being  described  in  detail  with  a  two 
half-circle  plot,  is  represented  with  sufficient  precision  with 
a  single  half-circle  plot  at  small  values  of  cathode  over 
potential;  this  corresponds  to  an  equivalent  network  for  this 
impedance  as  shown  in  Fig.  12b. 

In  particular,  the  following  tests  can  be  carried  out: 


1 .  Tests  at  continuously  increasing  frequency  (up  to  a  few 
kHz)  and  a  known  cell  temperature  T  and  gas  pressures 
and  small  values  of  cathode  overpotential  (low  currents); 
this  allows  to  determine  the  value  of  R„hm  as  the  value  of 
the  complex  impedance  at  the  left  intercept  of  the  half¬ 
circle  with  the  real  axis;  it  allows  also  the  determination 
of  Cdi  and  Rct:  the  experimental  behaviour  of  the 
complex  impedance  Z  is  to  be  replicated  by  the 
analytical  equation: 


Z(ico)  =  Rohm 


Rct 

1  +  i27t/RctCdi  ’ 


where  i  is  the  imaginary  unity. 

2.  Tests  at  different  currents  to  determine  the  behaviour  of 
/?ohm  with  the  current;  the  value  of  Rohm  is  determined  as 
the  left-intercept  of  measured  impedance  with  the  real 
axis. 

3.  Tests  as  those  described  in  the  previous  two  points,  but 
at  different  temperatures  and  gas  pressures  to  determine 
the  dependence  of  fi0hm  and  Cdi  on  these  two  quantities. 
In  this  work,  f?ohm  showed  linear  variation  with 
temperature  (Fig.  9),  while  the  parameter  Cdi  has  not 
shown  a  monotone  behaviour  with  the  temperature  and 
the  gas  pressures. 


3.3.  Parameters  oc#+  and  Th+ 

The  parameters  aH  ■  and  th+  appear  only  in  Eq.  (21).  If  the 
cell  is  subjected  to  a  sudden  reduction  in  the  cell  current,  the 
voltage  behaviour  is  of  the  type  shown  in  Fig.  13.  In  the 
figure,  the  part  of  the  transient  between  t\  and  t*  is  dominated 
by  the  inner  phenomena  related  to  gas  diffusion  and  double¬ 
layer  space  charge  modifications. 
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V  time  1 

If  the  cell  is  subjected  to  sudden  current  reductions  of 
constant  value  Aj  and  different  final  value  jf,  a  table  can 
be  built  in  which  values  of  exp(Avfc/fc)are  associated  to  the 
corresponding  values  of  jfk.  The  parameter  aH+  can  be 
identified  by  interpolating  the  values  exp(A vdb)  by  means 
of  the  following  function: 


/O'f)  = 


i  +  +  A /) 

1  +  aH+Jf 


(30) 


Considering  the  Eq.  (21),  in  the  interval  (f*,  ff)  the  voltage 
behaviour  is  an  exponential  function  having  as  time  constant 
th+  .  This  time  constant  can  therefore  be  inferred  from  the 
voltage  behaviour  as  that  shown  in  Fig.  13. 


3.4.  Parameters  %Ld,  r.gLd,  airef,  kh  y.h  ct2,  a2,rej;  k2, 

«3,  jrefi  ka,  rJ-o 


Between  t*  and  tf,  on  the  contrary,  the  phenomenon 
interpreted  in  the  model  as  proton  concentration  variations 
is  the  major  one.  If  this  is  assumed,  in  turn  b  and  p\  can  be 
assumed  constant  and  it  can  be  stated  jT  =  jf,  therefore, 


Jr  =  Jf  =J<A 


fP$P] 

bio[H+] 


=  expf-^  =  constant  =  K. 

[H+]0 


If  steady-state  voltage-current  density  curves  are  mea¬ 
sured  with  different  reactants  pressures  and  temperatures,  by 
means  of  an  error  minimisation  function  applied  to  Eq.  (27), 
a  set  of  numerical  values  can  be  found  for  each  of  the 
parameters  jol,  Am,  a\,  a2  and  v. 

Utilising  these  sets,  in  turn,  by  means  of  an  error  mini¬ 
misation  function  applied  to  the  equations  reported  in  the 
left  part  of  the  following  table,  the  numerical  values  of 
parameters  reported  in  its  right  part  can  be  determined: 


Equation 

Parameters  to  be  determined 

*  .  RT  r  psat  i  l 

(^d)2 

(rLd)2/egLd 

4coF  4 F  l(p  -  Psm)D]2  +  O13J 

egTd 

aiATX2  (h  h\ 

«l  =  «l.refP+Nai-r  exp  F —  -  —  j 

a  |  >ref,  kf,  OL  i,  «2 

Ps at,7ot3  fkl  k2  \ 

«2  =  «M-JV£exp(i--— j 

«2,ref>  k2,  Ct-3 

MA„=ja,a  p(-^  +  A;)pS 

Jref,  k„,  a0 

Then, 


[h+t 

[H+]0 

[H+r 

[H+]f 


m 

"  [H+], 


exp(!)> 

=  exp(^ 


Finally,  neglecting  the  proton  concentration  variation  in  the 
interval  (t-u  t*),  it  is 


[H+r 

[H+]f 


=  exp 


[H+^l+aH+J? 
[H+]f  1  +  %+/f 


1  +  (XH+J'f 


(29) 


Note  that  using  just  the  steady-state  cell  behaviour  it  is  not 
possible  to  separate  (t Ld)2/egLd  into  its  components  tLd  and 
egLd;  this  can  be  done  instead  by  resorting  to  the  dynamic 
cell  behaviour,  e.g.  the  response  to  a  step  current,  by  means 
of  a  trial  and  error  procedure. 


4.  Comparison  of  simulated  and  measurement  results 

4.1.  Experimental  set-up  and  model  parameter  values 

The  cell  under  test,  having  a  section  of  50  cm2,  has  been 
assembled  at  the  ENEA’s  labs.  The  cathode  had  a  catalyst 
E-TEK  20%  Pt/C  with  0.34  mg/cm2  of  Pt,  while  the  anode 
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Table  4 

Comparison  between  the  values  of  some  model  parameters  and  corresponding  values  reported  in  literature 


r Ld  (cm) 

(cm) 

Qi 

(mF/cm2) 

(mA/cm") 

U 

(cm) 

T 

% 

T  [°Cj 

Cathode 

Pt  loading 
[mg/cm2] 

Ref. 

0.1490 

0.0263 

41.71 

3.88 

- 

- 

- 

70 

0.34 

This  work 

- 

0.0104 

- 

10~5 

0.0260 

- 

0.40 

80 

- 

[6] 

- 

0.0066 

- 

90.00 

0.0300 

- 

0.22 

80 

0.12 

[7] 

- 

- 

- 

10.80 

0.0250 

- 

- 

50-70 

0.4 

[9] 

0.0850 

0.0020 

- 

5.9  x  104 

0.0170 

5 

0.12 

80 

0.10 

[11] 

0.2100 

0.0120 

13.70 

32.80 

0.0300 

7 

0.40 

80 

0.15 

[16] 

had  a  catalyst  realised  in  laboratory  with  20%  Pt/C  and 
0.29  mg/cm2  of  Pt.  The  diffusion  layer  had  300  |4m  of 
thickness  and  the  membrane  was  a  Nation  115  (125  |4m). 

The  electrodes-membrane  assembly  has  been  pressed 
between  two  graphite  plates  provided  with  serpentine  chan¬ 
nels.  Reactant  gases  were  humidified  Air  and  hydrogen. 

In  this  section  a  comparison  of  simulated  and  measure¬ 
ment  results  is  presented,  both  in  steady-state  and  dynamic 
conditions. 

The  simulations  are  carried  out  with  the  presented  model 
and  the  numerical  values  of  the  different  parameters  sum¬ 
marised  in  Table  3. 

A  comparison  of  values  of  all  parameters  in  Table  3  to 
values  of  the  same  parameters  which  can  be  found  in  literature 
has  not  been  possible;  in  fact  the  values  found  in  literature  has 
been  only  ones  of  the  followings:  r Ld,  egLd,  Cd|,  jref  and  Rmf. 

In  Section  2.6,  a  comparison  has  been  done  about  the 
membrane  conductivity  and  indirectly  about  the  parameter 

*ref. 

In  Table  4,  a  comparison  about  the  other  parameters  is 
reported.  As  it  can  be  seen,  the  values  of  parameters  iLd,  egLd, 
Cdl,  agree  to  ones  reported  in  literature,  while  nothing  can  be 
said  about  the  parameter  jTef,  considering  that  its  values 
reported  in  literature  are  in  heavy  disagreement  among  them. 


4.2.  Steady-state  cell  behaviour 

The  comparison  is  made  both  in  steady  state  and  dynamic 
conditions. 

In  Fig.  14,  a  comparison  of  steady-state  experimental  and 
simulated  cell  voltage  under  different  operating  conditions  is 
reported. 

The  figure  shows  a  good  agreement  of  data,  especially  in 
the  region  of  low  densities  (before  the  knee);  this  region  is 
the  one  having  the  maximum  interest  for  practical  applica¬ 
tions;  a  sufficient  agreement  is  however  present  also  even 
after  the  curve  knees. 

4.3.  Dynamic  cell  behaviour 

In  Figs.  15  and  16,  the  dynamic  response  to  current  steps 
of  the  model  is  compared  with  the  one  experimentally 
measured,  in  different  cell  operating  conditions. 

In  each  case  a  good  agreement  is  observed  on  the  whole 
transient. 

In  particular  in  the  case  of  Fig.  16,  showing  the  response 
to  an  instantaneous  interruption  of  the  cell  current,  the 
agreement  between  experimental  data  and  simulation  is 
excellent. 


Fig.  14.  Comparison  of  simulated  and  measured  cell  voltage  at  a  cell  temperature  of  70  °C  and  different  values  of  reactant  pressures  and  flows. 
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time  [s]  time  [s] 


Fig.  15.  Comparison  of  experimental  and  simulated  cell  dynamic  response,  considering  two  different  current  steps.  Cell  temperature  60  °C;  reactant 
pressures  3  bar;  air  and  H2  flows  1000  and  500  scc/min  respectively. 


5.  Conclusion 


continuity  equations : 


A  dynamical  model  of  proton  exchange  membrane  fuel 
cell,  based  on  physical  and  chemical  equations,  has  been 
proposed  and  described  in  detail. 

A  procedure  for  determining  the  numerical  value  of  the  22 
parameters  for  the  proposed  model  has  been  proposed,  the 
procedure  has  been  applied  to  a  real  cell,  and  the  values  have 
been  determined.  A  partial  comparison  of  these  values  has 
been  performed  with  similar  results  found  in  literature. 

The  behaviour  of  the  proposed  model,  used  with  values  of 
the  numerical  parameters  found  with  the  proposed  identifi¬ 
cation  procedure,  has  shown  a  good  agreement  with  experi¬ 
mental  results  both  in  static  and  dynamic  conditions. 

This  model  is  structured  so  that  it  can  be  easily  extended 
to  a  whole  stack  of  fuel  cells. 

It  is  believed  that  simplified  versions  of  this  model  are  to 
be  derived,  so  that  to  avoid  the  presence  of  partial  differ¬ 
ential  equations,  and  to  reduce  the  number  of  parameters  to 
be  identified.  Work  is  under  way  in  this  direction. 


Appendix  A.  Simplified  equations  for  gas 
diffusion  in  cathode 

The  gas  transport  in  cathode  is  described  by  the  following 
equations  considering  air  as  an  ideal  gas  (three-component 
gas): 
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Table  5 

Critical  properties  of  oxygen,  water  vapour  and  nitrogen 

Subscript  Species  M  Tc  (K)  pc  (bar) 

1  02  31.999  154.35  50.40 

2  H20  18.015  647.30  221.2 

3  N2  28.013  126.05  33.93 


the  subscripts  1,  2,  and  3  in  previous  equations  are  rela¬ 
ted,  respectively,  to  components  oxygen,  water  vapour  and 
nitrogen. 

The  pressure-diffusivity  products  (barcm2/s)  are  esti¬ 
mated  [13]  from  critical  temperature,  critical  pressure, 
and  mass  of  the  respective  components  (Table  5): 


=  D12(T)  =  3.55x10-4 
=  4.281  x  10“7r2-334 


/j_noW/3  r 

\Mi  M2){TciTc2f1^: 


p2>n  =  Dn(T)  =  1.882  xl(T3 
=  5.338  x  10_57’L5/exp| 


Mi  M3 

'83.63' 


■)  /(4.Gi3  )T3/2 


pS>23  =  D23  (T)  =  3.687  x  10“4 
=  4.477  x  l0-7r2-334 


Vl  n~faw)1/3 

\m2  m3)  (rc2rc3)0750 


j-2.334 


(A.  11) 


A  comparison  of  (A.9)  and  (A.l  1)  yields  D23  =  D\2  (D23l 
D\2  =  1.046);  for  sake  of  simplicity,  in  this  work  the 
following  is  assumed: 

£>23  =  £>12  (A.  12) 

As  reported  in  Section  2.3  the  water  vapour  contained  in  the 
gas  reactants  is  considered  in  equilibrium  with  the  surround¬ 
ing  liquid  phase,  and  its  saturation  pressure  psat  used  here 
was  fitted  to  the  following  expression  according  to  tabulated 
values  in  [19]: 

Psat(£)(bar)  =  Aexp^-^j  =7.714  x  105  exp^-50^9^ . 

(A.  13) 


Using  the  assumptions  3  and  4  in  Section  2.3  and  the 
Eqs.  (A.7),  (A.8)  and  (A.  12)  yields  relatively  at  Eqs.  (A.2) 
and  (A.  5) 


dN2_ 


(A.  14) 


N2  =  N2(t)  =  [Ni  (y,  t)  +  N3(y,  f)]  (A.  15) 

Pc  -  P2 

At  diffusion  layer/catalyst  layer  interface  (Fig.  1,  y  =  Ld)  it 
is  assumed  that  the  instantaneous  oxygen  flux  is  equal  to 
total  oxygen  flux  consumed  in  catalyst  layer: 


Ai(£d,t) 


(A.  16) 


the  oxygen  storage  in  catalyst  layer  is  considered  negligible. 


As  nitrogen  is  inert  and  its  solubility  in  water  is  very  little 
it  is  assumed  that  the  Nitrogen  flux  in  y  =  £d  is  null: 

A3(£d,f)=0.  (A.  17) 

Then  using  Eqs.  (A.16)  and  (A.17)  the  Eq.  (A.15)  can  be 
written  as 


N2  =  N2(t)  =  — V—  £=.  (A.  18) 

Pc-Pi^F 

Substituting  the  Eq.  (A.  18)  in  Eq.  (A.4),  deriving  the  result¬ 
ing  equation  with  respect  to  y  and  combining  the  result  with 
Eq.  (A.l)  the  following  single  partial  differential  equation  is 
obtained 


in  which  ^  is  the  dimensionless  distance  y/£d  and 


0>  T2£2(p2/£>12  +  {Pc  ~  P2)/£>13) 

1 

*2£dOsat/£>12  +  (Pc-Pm)/D13) 

(A.  20) 

RT  RT 

eg£d(pc  -  Pi)  eg£d(/>c  -  Pm) ' 

(A.21) 

To  solve  the  partial  differential  Eq.  (A.  19)  it  is  necessary  to 
define  two  boundary  conditions.  At  flow  channels/diffusion 
layer  interface  (Fig.  1,  y  =  0)  the  condition 

pi  (0,  t)  =  pio  =  ^  =  constant  (A.22) 

1  +P 


applies,  assuming  constant  oxygen  partial  pressure  in  flow 
channels.  The  parameter  [1  is  the  nitrogen-oxygen  mole  ratio 
of  inlet  air. 

The  second  boundary  condition  is  obtained  evaluating  the 
Eq.  (A.4)  in  ^  =  1  and  using  the  conditions  (A.16),  (A.17) 
and  (A.  18): 

^(ht)  =  ^\pi(l,t)+p^-pc}  (A.23) 
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